Introduction

Here is my analysis of data/Carseats_org.csv.

Configuration

library(leaps)
library(stringr)
library(caret)
library(ggplot2)
library(DataExplorer)
library(dplyr)
library(ggExtra)
library(RColorBrewer)
library(plotly)
library(corrplot)
library(htmltools)
library(MASS)

System Information

Due to the large number of libraries in use I have provided system information.

sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /anaconda3/lib/R/lib/libRblas.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] MASS_7.3-51.1      htmltools_0.3.6    corrplot_0.84      plotly_4.8.0       RColorBrewer_1.1-2 ggExtra_0.8       
 [7] dplyr_0.8.0.1      DataExplorer_0.7.0 caret_6.0-81       ggplot2_3.1.0      lattice_0.20-38    stringr_1.4.0     
[13] leaps_3.0         

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0         lubridate_1.7.4    tidyr_0.8.2        class_7.3-15       assertthat_0.2.0   digest_0.6.18     
 [7] ipred_0.9-8        foreach_1.4.4      mime_0.6           R6_2.4.0           plyr_1.8.4         stats4_3.5.1      
[13] evaluate_0.12      httr_1.4.0         pillar_1.3.1       rlang_0.3.1        lazyeval_0.2.1     rstudioapi_0.9.0  
[19] data.table_1.12.0  miniUI_0.1.1.1     rpart_4.1-13       Matrix_1.2-15      rmarkdown_1.11     labeling_0.3      
[25] splines_3.5.1      gower_0.1.2        htmlwidgets_1.3    igraph_1.2.4       munsell_0.5.0      shiny_1.2.0       
[31] compiler_3.5.1     httpuv_1.4.5.1     xfun_0.4           pkgconfig_2.0.2    base64enc_0.1-3    nnet_7.3-12       
[37] tidyselect_0.2.5   tibble_2.0.1       gridExtra_2.3      prodlim_2018.04.18 codetools_0.2-16   viridisLite_0.3.0 
[43] later_0.8.0        crayon_1.3.4       withr_2.1.2        recipes_0.1.4      ModelMetrics_1.1.0 grid_3.5.1        
[49] xtable_1.8-3       nlme_3.1-137       jsonlite_1.6       gtable_0.2.0       magrittr_1.5       scales_1.0.0      
[55] stringi_1.3.1      reshape2_1.4.3     promises_1.0.1     timeDate_3043.102  generics_0.0.2     lava_1.6.5        
[61] iterators_1.0.10   tools_3.5.1        glue_1.3.0         purrr_0.2.5        crosstalk_1.0.0    networkD3_0.4     
[67] rsconnect_0.8.13   parallel_3.5.1     survival_2.43-3    yaml_2.2.0         colorspace_1.4-0   knitr_1.21        
sapply(c('repr', 'IRdisplay', 'IRkernel'), function(p) paste(packageVersion(p)))
     repr IRdisplay  IRkernel 
 "0.19.2"   "0.7.0"  "0.8.15" 

Data

I load the data into r, and drop the “ID” column.

carseats <- read.csv("data/Carseats_org.csv", header = T, stringsAsFactors = T)
drops <- c("X")
carseats <- carseats[ , !(names(carseats) %in% drops)]
head(carseats, 10)

Here I create two new data frames to manage numeric and categorical data.

# get vectors of continuous and categorical cols
nums <- dplyr::select_if(carseats, is.numeric)
cats <- dplyr::select_if(carseats, is.factor)
nums[sample(nrow(nums), 10), ]
cats[sample(nrow(cats), 10), ]

Let’s get some quick summaries of each:

print('Numeric Summaries')
[1] "Numeric Summaries"
summary(nums)
     Sales          CompPrice       Income        Advertising       Population        Price            Age       
 Min.   : 0.000   Min.   : 77   Min.   : 21.00   Min.   : 0.000   Min.   : 10.0   Min.   : 24.0   Min.   :25.00  
 1st Qu.: 5.390   1st Qu.:115   1st Qu.: 42.75   1st Qu.: 0.000   1st Qu.:139.0   1st Qu.:100.0   1st Qu.:39.75  
 Median : 7.490   Median :125   Median : 69.00   Median : 5.000   Median :272.0   Median :117.0   Median :54.50  
 Mean   : 7.496   Mean   :125   Mean   : 68.66   Mean   : 6.635   Mean   :264.8   Mean   :115.8   Mean   :53.32  
 3rd Qu.: 9.320   3rd Qu.:135   3rd Qu.: 91.00   3rd Qu.:12.000   3rd Qu.:398.5   3rd Qu.:131.0   3rd Qu.:66.00  
 Max.   :16.270   Max.   :175   Max.   :120.00   Max.   :29.000   Max.   :509.0   Max.   :191.0   Max.   :80.00  
   Education   
 Min.   :10.0  
 1st Qu.:12.0  
 Median :14.0  
 Mean   :13.9  
 3rd Qu.:16.0  
 Max.   :18.0  
print('Categorical Summaries')
[1] "Categorical Summaries"
summary(cats)
  ShelveLoc   Urban       US     
 Bad   : 96   No :118   No :142  
 Good  : 85   Yes:282   Yes:258  
 Medium:219                      
str(nums)
'data.frame':   400 obs. of  8 variables:
 $ Sales      : num  9.5 11.22 10.06 7.4 4.15 ...
 $ CompPrice  : int  138 111 113 117 141 124 115 136 132 132 ...
 $ Income     : int  73 48 35 100 64 113 105 81 110 113 ...
 $ Advertising: int  11 16 10 4 3 13 0 15 0 0 ...
 $ Population : int  276 260 269 466 340 501 45 425 108 131 ...
 $ Price      : int  120 83 80 97 128 72 108 120 124 124 ...
 $ Age        : int  42 65 59 55 38 78 71 67 76 76 ...
 $ Education  : int  17 10 12 14 13 16 15 10 10 17 ...
str(cats)
'data.frame':   400 obs. of  3 variables:
 $ ShelveLoc: Factor w/ 3 levels "Bad","Good","Medium": 1 2 3 3 1 1 3 2 3 3 ...
 $ Urban    : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
 $ US       : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...

Data Dimensionality

This command is to inspect the different data types in the data.

str(carseats)
'data.frame':   400 obs. of  11 variables:
 $ Sales      : num  9.5 11.22 10.06 7.4 4.15 ...
 $ CompPrice  : int  138 111 113 117 141 124 115 136 132 132 ...
 $ Income     : int  73 48 35 100 64 113 105 81 110 113 ...
 $ Advertising: int  11 16 10 4 3 13 0 15 0 0 ...
 $ Population : int  276 260 269 466 340 501 45 425 108 131 ...
 $ Price      : int  120 83 80 97 128 72 108 120 124 124 ...
 $ ShelveLoc  : Factor w/ 3 levels "Bad","Good","Medium": 1 2 3 3 1 1 3 2 3 3 ...
 $ Age        : int  42 65 59 55 38 78 71 67 76 76 ...
 $ Education  : int  17 10 12 14 13 16 15 10 10 17 ...
 $ Urban      : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
 $ US         : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...
print("")
[1] ""
print(paste('Number of Columns:', ncol(carseats)))
[1] "Number of Columns: 11"
print(paste('Number of Numeric Columns:', ncol(nums)))
[1] "Number of Numeric Columns: 8"
print(paste('Number of Categorical Columns:', ncol(cats)))
[1] "Number of Categorical Columns: 3"
dim(carseats)
[1] 400  11
dim(nums)
[1] 400   8
dim(cats)
[1] 400   3

Here’s a quick way to examine general properties of the data:

DataExplorer::introduce(data=carseats)

Finally, I want to look at the first and last rows of the data set. Just to be safe:

head(carseats, 2)
tail(carseats, 2)
head(nums, 2)
tail(nums, 2)
head(cats, 2)
tail(cats, 2)

Numeric Plotting

I start out with a few general scatter plots.

plot_ly(data=carseats, 
        x=~Age, 
        y=~Sales, 
        mode = 'markers', 
        type = 'scatter', 
        color=~ShelveLoc) %>%
            layout(title = "Age, Shelf Location, and Sales Scatter Plot", width=900) 
Specifying width/height in layout() is now deprecated.
Please specify in ggplotly() or plot_ly()

This plot below shows good separation and a weak linear trend. These variables are worth investigating.

plot_ly(data=carseats, 
        x=~Price, 
        y=~Sales, 
        mode = 'markers', 
        type = 'scatter', 
        color=~ShelveLoc) %>%
            layout(title = "Price, Shelf Location, and Sales Scatter Plot", width=900)
Specifying width/height in layout() is now deprecated.
Please specify in ggplotly() or plot_ly()

Here we inspect the density of the ‘Price vs Sales’ relationship:

plot_ly(data=carseats, 
        x=~Price, 
        y=~Sales, 
        mode = 'markers', 
        size = ~Price,
        type = 'scatter',
        colors = "Dark2",
        alpha = .6) %>%
            layout(title = "Price, US, and Sales Scatter Plot", width=900) 
Specifying width/height in layout() is now deprecated.
Please specify in ggplotly() or plot_ly()`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.

Normalization

I choose to normalize the numeric data in order to be able to plot each variable on the same scale. This will allow me to investigate the variation of each predictor relative to Sales.

preObj <- preProcess(nums, method=c("center", "scale"))
scaled.nums <- predict(preObj, nums)
head(scaled.nums, 2)
tail(scaled.nums, 2)
str(scaled.nums)
'data.frame':   400 obs. of  8 variables:
 $ Sales      : num  0.7095 1.3185 0.9078 -0.0341 -1.1849 ...
 $ CompPrice  : num  0.849 -0.911 -0.781 -0.52 1.045 ...
 $ Income     : num  0.155 -0.738 -1.203 1.12 -0.166 ...
 $ Advertising: num  0.656 1.408 0.506 -0.396 -0.547 ...
 $ Population : num  0.0757 -0.0328 0.0282 1.3649 0.51 ...
 $ Price      : num  0.178 -1.385 -1.512 -0.794 0.515 ...
 $ Age        : num  -0.699 0.721 0.35 0.104 -0.946 ...
 $ Education  : num  1.183 -1.4882 -0.725 0.0382 -0.3434 ...
print("")
[1] ""
summary(scaled.nums)
     Sales            CompPrice            Income          Advertising        Population           Price         
 Min.   :-2.65440   Min.   :-3.12856   Min.   :-1.70290   Min.   :-0.9977   Min.   :-1.72918   Min.   :-3.87702  
 1st Qu.:-0.74584   1st Qu.:-0.65049   1st Qu.:-0.92573   1st Qu.:-0.9977   1st Qu.:-0.85387   1st Qu.:-0.66711  
 Median :-0.00224   Median : 0.00163   Median : 0.01224   Median :-0.2459   Median : 0.04858   Median : 0.05089  
 Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.00000  
 3rd Qu.: 0.64575   3rd Qu.: 0.65375   3rd Qu.: 0.79834   3rd Qu.: 0.8067   3rd Qu.: 0.90693   3rd Qu.: 0.64219  
 Max.   : 3.10670   Max.   : 3.26225   Max.   : 1.83458   Max.   : 3.3630   Max.   : 1.65671   Max.   : 3.17633  
      Age             Education       
 Min.   :-1.74827   Min.   :-1.48825  
 1st Qu.:-0.83779   1st Qu.:-0.72504  
 Median : 0.07268   Median : 0.03816  
 Mean   : 0.00000   Mean   : 0.00000  
 3rd Qu.: 0.78255   3rd Qu.: 0.80137  
 Max.   : 1.64673   Max.   : 1.56457  

Distributions

Here are scaled distributions (histograms and density plots) for each numeric variable, including Sales. The variables relating to money ($) tend to be approximately normal. Many other variables tend to be approximately uniform, which does not bode well for their predictive power.

scaled.nums %>%
    tidyr::gather() %>% 
        ggplot(aes(x=value,y=..density..))+
            ggtitle('Distributions of Continous Variables (scaled)') +
            facet_wrap(~ key, scales = "free") +
            geom_histogram(fill=I("orange"), col=I("black"), bins = 50) +
            facet_wrap(~ key, scales = "free") +
            geom_density(color="blue", fill='light blue', alpha = 0.4)

Here we plot all numeric variables against their distributions. This is just another way to examine the information shown above.

scaled.nums %>%
    tidyr::gather() %>% 
        plot_ly(x=~key, y=~value,
                type = "box", 
                boxpoints = "all", 
                jitter = 0.4,
                pointpos = 0,
                color = ~key,
                colors = "Dark2") %>%
                      subplot(shareX = TRUE)  %>%
                            layout(title = "Numeric Variable Distributions (scaled)", 
                                  yaxis=list(title='Standard Deviation'),
                                  xaxis=list(title='Variable'),
                                  autosize=FALSE,
                                  width=900,
                                  height=500)
Specifying width/height in layout() is now deprecated.
Please specify in ggplotly() or plot_ly()

Scatterplots

Here we plot all numeric variables against Sales (scaled). This allows us to investigate possible linear relationships between that variable and Sales. As shown below, only ‘Price’ appears to have a linear relationship worth investigating. This took me so long to figure out.

numeric.scatterplots <- htmltools::tagList()
count = 1
for (i in names(scaled.nums[,-1])) {
  
  numeric.scatterplots[[count]] <- plot_ly(scaled.nums, x=scaled.nums[,i], y=scaled.nums$Sales,
                    colors = 'RdYlGn',
                    mode = 'markers',
                    type = 'scatter',
                    size = scaled.nums$Sales^2,
                    color = scaled.nums$Sales,
                    marker = list(line = list(color = 'black',width = 2)),
                    name=paste(i)) %>%
          layout(title = paste(i, "vs Sales (scaled)<br>Size=Sales^2"),
                               yaxis=list(title='Sales'),
                               xaxis=list(title=i),
                               showlegend = FALSE)
  
  count = count + 1
  
}
numeric.scatterplots
`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.

By adding naive regression lines to a few scatter plots we can confirm our suspicions:

fit.Pop <- lm(Sales ~ Population, data = scaled.nums)
fit.Age <- lm(Sales ~ Age, data = scaled.nums)
fit.CompPrice <- lm(Sales ~ CompPrice, data = scaled.nums)
fit.Price <- lm(Sales ~ Price, data = scaled.nums)
regression.scatterplots <- htmltools::tagList()
regression.scatterplots[[1]] <-  plot_ly(scaled.nums, 
          x = ~Population,
          name = 'Population vs Sales Regression Line') %>% 
              add_markers(y = ~Sales,
                    name = 'Population vs Sales Observations') %>% 
                            add_lines(x = ~Population, 
                                  y = fitted(fit.Pop)) %>%
                                        layout(title = "Population vs Sales",
                                               yaxis=list(title='Sales',
                                               xaxis=list(title='Population')),
                                               showlegend = FALSE)
                                              
           
regression.scatterplots[[2]] <-  plot_ly(scaled.nums, 
          x = ~Age,
          name = 'Age vs Sales Regression Line') %>% 
              add_markers(y = ~Sales,
                    name = 'Age vs Sales Observations') %>% 
                            add_lines(x = ~Age, 
                                  y = fitted(fit.Age)) %>%
                                        layout(title = "Age vs Sales",
                                               yaxis=list(title='Sales',
                                               xaxis=list(title='Age')),
                                               showlegend = FALSE)
regression.scatterplots[[3]] <-  plot_ly(scaled.nums, 
          x = ~CompPrice,
          name = 'CompPrice vs Sales Regression Line') %>% 
              add_markers(y = ~Sales,
                    name = 'CompPrice vs Sales Observations') %>% 
                            add_lines(x = ~CompPrice, 
                                  y = fitted(fit.CompPrice)) %>%
                                        layout(title = "CompPrice vs Sales",
                                               yaxis=list(title='Sales',
                                               xaxis=list(title='CompPrice')),
                                               showlegend = FALSE)
regression.scatterplots[[4]] <-  plot_ly(scaled.nums, 
          x = ~Price,
          name = 'Price vs Sales Regression Line') %>% 
              add_markers(y = ~Sales,
                    name = 'Price vs Sales Observations') %>% 
                            add_lines(x = ~Price, 
                                  y = fitted(fit.Price)) %>%
                                        layout(title = "Price vs Sales",
                                               yaxis=list(title='Sales',
                                               xaxis=list(title='Price')),
                                               showlegend = FALSE)
regression.scatterplots

Let’s compare the slopes:

scaled.nums %>%
    plot_ly(y = ~Sales) %>%
      add_lines(x= ~Population, y = fitted(fit.Pop), 
                name = "fit.Pop slope", line = list(shape = "linear")) %>%
      add_lines(x= ~Age, y = fitted(fit.Age), 
                name = "fit.Age slope", line = list(shape = "linear")) %>%
      add_lines(x= ~CompPrice, y = fitted(fit.CompPrice), 
                name = "fit.CompPrice slope", line = list(shape = "linear")) %>%
      add_lines(x= ~Price, y = fitted(fit.Price), 
                name = "fit.Price slope", line = list(shape = "linear")) %>%
            
    layout(title = "Regression Lines vs Sales (scaled)",
           autosize=FALSE,
           width=900,
           yaxis=list(title='Sales'),
           xaxis=list(title='Scaled Numeric Variable'))
Specifying width/height in layout() is now deprecated.
Please specify in ggplotly() or plot_ly()

Here’s a pretty graphic that doesn’t help me understand anything about the data.

y = scaled.nums$Sales
x = scaled.nums$Price
s <- subplot(
  plot_ly(x = x, color = I("black"), type = 'histogram'), 
  plotly_empty(), 
  plot_ly(x = x, y = y, type = 'histogram2dcontour', showscale = F), 
  plot_ly(y = y, color = I("black"), type = 'histogram'),
  nrows = 2, heights = c(0.2, 0.8), widths = c(0.8, 0.2),
  shareX = TRUE, shareY = TRUE, titleX = FALSE, titleY = FALSE)
No trace type specified and no positional attributes specifiedNo trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
layout(s, showlegend = FALSE, autosize=FALSE,
           width=700,
           height=500, 
           yaxis=list(title='Sales'),
           xaxis=list(title='Price'))
Specifying width/height in layout() is now deprecated.
Please specify in ggplotly() or plot_ly()

Categorical Plotting

First, let’s create a data frame that we can use:

categorical.by.sales = cbind(Sales = scaled.nums$Sales, cats)
str(categorical.by.sales)
'data.frame':   400 obs. of  4 variables:
 $ Sales    : num  0.7095 1.3185 0.9078 -0.0341 -1.1849 ...
 $ ShelveLoc: Factor w/ 3 levels "Bad","Good","Medium": 1 2 3 3 1 1 3 2 3 3 ...
 $ Urban    : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
 $ US       : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...

Here we can see all categorical by Sales. We suspected that ‘ShelveLoc’ would be important based on one of the early scatter plots. It seems that this is the case.

categorical.boxplots <- htmltools::tagList()
count = 1
for (i in names(categorical.by.sales[,-1])) {
  
  categorical.boxplots[[count]] <- plot_ly(categorical.by.sales, x=categorical.by.sales[,i], y=categorical.by.sales$Sales,
                        type = "box", 
                        boxpoints = "all",
                        jitter = .2,
                        pointpos = 0,
                        color =categorical.by.sales[,i],
                        colors='Set1', 
                        name=paste(i)) %>%
                            layout(title = paste(i, "vs Sales (scaled)"),
                                  showlegend = TRUE,
                                  yaxis=list(title='Sales Standard Deviation'),
                                  xaxis=list(title=i))
  
  count = count + 1
  
}
categorical.boxplots

Here’s the same thing, but more musically:

categorical.violins <- htmltools::tagList()
count = 1
for (i in names(categorical.by.sales[,-1])) {
  
  categorical.violins[[count]] <- plot_ly(categorical.by.sales, x=categorical.by.sales[,i], y=categorical.by.sales$Sales,
                        split = categorical.by.sales[,i],
                        type = 'violin',
                        colors='Set1', 
                        name=paste(i),
                        box = list(visible = TRUE),
                        meanline = list(visible = TRUE)) %>% 
                            layout(xaxis = list(title = "US"),
                                yaxis = list(title = "Sales",zeroline = FALSE))
  
  count = count + 1
  
}
categorical.violins

Linear Regression

First, let’s merge the data set into a single data frame

scaled.merged <- cbind(categorical.by.sales[,-1], scaled.nums)
str(scaled.merged)
'data.frame':   400 obs. of  11 variables:
 $ ShelveLoc  : Factor w/ 3 levels "Bad","Good","Medium": 1 2 3 3 1 1 3 2 3 3 ...
 $ Urban      : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
 $ US         : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...
 $ Sales      : num  0.7095 1.3185 0.9078 -0.0341 -1.1849 ...
 $ CompPrice  : num  0.849 -0.911 -0.781 -0.52 1.045 ...
 $ Income     : num  0.155 -0.738 -1.203 1.12 -0.166 ...
 $ Advertising: num  0.656 1.408 0.506 -0.396 -0.547 ...
 $ Population : num  0.0757 -0.0328 0.0282 1.3649 0.51 ...
 $ Price      : num  0.178 -1.385 -1.512 -0.794 0.515 ...
 $ Age        : num  -0.699 0.721 0.35 0.104 -0.946 ...
 $ Education  : num  1.183 -1.4882 -0.725 0.0382 -0.3434 ...
head(nums, 2)
tail(nums, 2)
head(scaled.merged, 2)
tail(scaled.merged, 2)

First, let’s look at some things that may give us trouble. Luckily it looks like the only serious correlation is with our dependent variable. We’ll want to watch the ‘Price’ vs ‘CompPrice’ relationship.

res <- cor(scaled.nums)
corrplot(res, type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45)

It appears that residuals are roughly symmetrical around 0. That’s strange. Mostly due to a relatively poor overall fit. Note how close to zero most of the coefficient estimates are.

simple.lm <- lm(Sales~., data=scaled.merged)
simple.summary <- summary(simple.lm)
print(simple.summary)

Call:
lm(formula = Sales ~ ., data = scaled.merged)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.01598 -0.24463  0.00748  0.23496  1.20797 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -0.73292    0.05999 -12.217  < 2e-16 ***
ShelveLocGood    1.71742    0.05422  31.678  < 2e-16 ***
ShelveLocMedium  0.69286    0.04465  15.516  < 2e-16 ***
UrbanYes         0.04351    0.04000   1.088    0.277    
USYes           -0.06519    0.05306  -1.229    0.220    
CompPrice        0.50397    0.02252  22.378  < 2e-16 ***
Income           0.15660    0.01828   8.565 2.58e-16 ***
Advertising      0.28987    0.02619  11.066  < 2e-16 ***
Population       0.01085    0.01933   0.561    0.575    
Price           -0.79946    0.02239 -35.700  < 2e-16 ***
Age             -0.26413    0.01825 -14.472  < 2e-16 ***
Education       -0.01958    0.01830  -1.070    0.285    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3608 on 388 degrees of freedom
Multiple R-squared:  0.8734,    Adjusted R-squared:  0.8698 
F-statistic: 243.4 on 11 and 388 DF,  p-value: < 2.2e-16

Linear Models and Subsets

Let’s do the same thing, but control the subsets using leaps

regfit.full=regsubsets(Sales~., data=scaled.merged, nvmax=5)
reg.summary=summary(regfit.full)
print(reg.summary)
Subset selection object
Call: regsubsets.formula(Sales ~ ., data = scaled.merged, nvmax = 5)
11 Variables  (and intercept)
                Forced in Forced out
ShelveLocGood       FALSE      FALSE
ShelveLocMedium     FALSE      FALSE
UrbanYes            FALSE      FALSE
USYes               FALSE      FALSE
CompPrice           FALSE      FALSE
Income              FALSE      FALSE
Advertising         FALSE      FALSE
Population          FALSE      FALSE
Price               FALSE      FALSE
Age                 FALSE      FALSE
Education           FALSE      FALSE
1 subsets of each size up to 5
Selection Algorithm: exhaustive
         ShelveLocGood ShelveLocMedium UrbanYes USYes CompPrice Income Advertising Population Price Age Education
1  ( 1 ) "*"           " "             " "      " "   " "       " "    " "         " "        " "   " " " "      
2  ( 1 ) "*"           " "             " "      " "   " "       " "    " "         " "        "*"   " " " "      
3  ( 1 ) "*"           " "             " "      " "   "*"       " "    " "         " "        "*"   " " " "      
4  ( 1 ) "*"           " "             " "      " "   "*"       " "    "*"         " "        "*"   " " " "      
5  ( 1 ) "*"           "*"             " "      " "   "*"       " "    "*"         " "        "*"   " " " "      

We’ll just take code straight from the example on Canvas…

par(mfrow=c(2,2))
plot(reg.summary$rss,xlab="Number of Variables",ylab="RSS",type="l")
plot(reg.summary$adjr2,xlab="Number of Variables",ylab="Adjusted RSq",type="l")
plot(reg.summary$cp,xlab="Number of Variables",ylab="Cp",type='l')
plot(reg.summary$bic,xlab="Number of Variables",ylab="BIC",type='l')

plot(regfit.full,scale="r2")

plot(regfit.full,scale="adjr2")

plot(regfit.full,scale="Cp")

plot(regfit.full,scale="bic")

Interaction Terms

Here we define a new model with some interaction terms: a. Income and Advertising b. Income and CompPrice c. Price and Age

interaction.lm <- lm(Sales~. + Income*Advertising + Income*CompPrice + Price*Age, data=scaled.merged)
interaction.summary <- summary(interaction.lm)
print(interaction.summary)

Call:
lm(formula = Sales ~ . + Income * Advertising + Income * CompPrice + 
    Price * Age, data = scaled.merged)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.04967 -0.23955 -0.00936  0.24591  1.19774 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.730051   0.059236 -12.324  < 2e-16 ***
ShelveLocGood       1.714324   0.053519  32.032  < 2e-16 ***
ShelveLocMedium     0.680563   0.044177  15.405  < 2e-16 ***
UrbanYes            0.045257   0.039380   1.149  0.25117    
USYes              -0.069622   0.052329  -1.330  0.18415    
CompPrice           0.507567   0.022131  22.935  < 2e-16 ***
Income              0.152164   0.018025   8.442 6.43e-16 ***
Advertising         0.289455   0.025756  11.238  < 2e-16 ***
Population          0.007224   0.018985   0.381  0.70377    
Price              -0.797242   0.022003 -36.233  < 2e-16 ***
Age                -0.260783   0.017978 -14.505  < 2e-16 ***
Education          -0.024248   0.018063  -1.342  0.18027    
Income:Advertising  0.046378   0.018170   2.552  0.01108 *  
CompPrice:Income   -0.059184   0.018911  -3.130  0.00188 ** 
Price:Age           0.013171   0.017912   0.735  0.46259    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3538 on 385 degrees of freedom
Multiple R-squared:  0.8792,    Adjusted R-squared:  0.8748 
F-statistic: 200.1 on 14 and 385 DF,  p-value: < 2.2e-16
interaction.lm.subsets <- regsubsets(Sales~. + Income*Advertising + Income*CompPrice + Price*Age, 
                                     data=scaled.merged, nvmax=5)
interaction.subsets.summary <- summary(interaction.lm.subsets)
print(interaction.subsets.summary)
Subset selection object
Call: regsubsets.formula(Sales ~ . + Income * Advertising + Income * 
    CompPrice + Price * Age, data = scaled.merged, nvmax = 5)
14 Variables  (and intercept)
                   Forced in Forced out
ShelveLocGood          FALSE      FALSE
ShelveLocMedium        FALSE      FALSE
UrbanYes               FALSE      FALSE
USYes                  FALSE      FALSE
CompPrice              FALSE      FALSE
Income                 FALSE      FALSE
Advertising            FALSE      FALSE
Population             FALSE      FALSE
Price                  FALSE      FALSE
Age                    FALSE      FALSE
Education              FALSE      FALSE
Income:Advertising     FALSE      FALSE
CompPrice:Income       FALSE      FALSE
Price:Age              FALSE      FALSE
1 subsets of each size up to 5
Selection Algorithm: exhaustive
         ShelveLocGood ShelveLocMedium UrbanYes USYes CompPrice Income Advertising Population Price Age Education
1  ( 1 ) "*"           " "             " "      " "   " "       " "    " "         " "        " "   " " " "      
2  ( 1 ) "*"           " "             " "      " "   " "       " "    " "         " "        "*"   " " " "      
3  ( 1 ) "*"           " "             " "      " "   "*"       " "    " "         " "        "*"   " " " "      
4  ( 1 ) "*"           " "             " "      " "   "*"       " "    "*"         " "        "*"   " " " "      
5  ( 1 ) "*"           "*"             " "      " "   "*"       " "    "*"         " "        "*"   " " " "      
         Income:Advertising CompPrice:Income Price:Age
1  ( 1 ) " "                " "              " "      
2  ( 1 ) " "                " "              " "      
3  ( 1 ) " "                " "              " "      
4  ( 1 ) " "                " "              " "      
5  ( 1 ) " "                " "              " "      

Variable Significance

Below we print the coefficients for the 5th model using the default model selection criteria. All coefficients are relatively small, as we would expect from the EDA above. This pretty much confirms what I would have guessed by looking at the data against sales. We still want to watch out for confounding between ‘Price’ and ‘CompPrice.’

coef(interaction.lm.subsets, 5)
    (Intercept)   ShelveLocGood ShelveLocMedium       CompPrice     Advertising           Price 
     -0.6995348       1.6746198       0.6277225       0.5090177       0.2832405      -0.7846476 

Second Interaction Model

First, drop columns unneeded from analysis:

scaled.merged.slim <- scaled.merged[ , -which(names(scaled.merged) %in% c("US","Urban"))]

A few hyper parameters we’d like to be consistent for all models

nvmax <- 3

Forward Selection:

interaction.subset.fwd <- regsubsets(Sales~. + Income*Advertising + Income*CompPrice + Income*Age,
                                     data=scaled.merged.slim, 
                                     nvmax = nvmax,
                                     method="forward")
fwd.subset.summary <- summary(interaction.subset.fwd)
coef(interaction.subset.fwd, 1:nvmax)
[[1]]
  (Intercept) ShelveLocGood 
    -0.259671      1.221981 

[[2]]
  (Intercept) ShelveLocGood         Price 
   -0.2708256     1.2744733    -0.4688868 

[[3]]
  (Intercept) ShelveLocGood     CompPrice         Price 
   -0.2709362     1.2749938     0.4932455    -0.7573701 

Backward Selection:

This is really strange. I can’t seem to find any documentation about this, but it appears that this model is actually ‘forward.’

interaction.subset.bk <- regsubsets(Sales~. + Income*Advertising + Income*CompPrice + Income*Age,
                                     data=scaled.merged.slim, 
                                     nvmax = nvmax,
                                     method="backward")
bk.subset.summary <- summary(interaction.subset.bk)
coef(interaction.subset.bk, 1:nvmax)
[[1]]
  (Intercept) ShelveLocGood 
    -0.259671      1.221981 

[[2]]
  (Intercept) ShelveLocGood         Price 
   -0.2708256     1.2744733    -0.4688868 

[[3]]
  (Intercept) ShelveLocGood     CompPrice         Price 
   -0.2709362     1.2749938     0.4932455    -0.7573701 

Exhasutive

interaction.subset.ex <- regsubsets(Sales~. + Income*Advertising + Income*CompPrice + Income*Age,
                                     data=scaled.merged.slim, 
                                     nvmax = nvmax,
                                     method="exhaustive")
ex.subset.summary <- summary(interaction.subset.ex)
coef(interaction.subset.ex, 1:nvmax)
[[1]]
  (Intercept) ShelveLocGood 
    -0.259671      1.221981 

[[2]]
  (Intercept) ShelveLocGood         Price 
   -0.2708256     1.2744733    -0.4688868 

[[3]]
  (Intercept) ShelveLocGood     CompPrice         Price 
   -0.2709362     1.2749938     0.4932455    -0.7573701 

This is a list that may come in handy.

model.list = list(list("Forward", interaction.subset.fwd, fwd.subset.summary), 
                  list("Backward", interaction.subset.bk, bk.subset.summary), 
                  list("Exhaustive", interaction.subset.ex, ex.subset.summary))

Evaluation Metric Plotting

It is interesting to see that each model selected the same variables, in the same order.

par(mfrow=c(2,2))
plot(fwd.subset.summary$rss,xlab="Number of Variables (forward)",ylab="RSS",type="l")
plot(fwd.subset.summary$adjr2,xlab="Number of Variables (forward)",ylab="Adjusted RSq",type="l")
plot(fwd.subset.summary$cp,xlab="Number of Variables (forward)",ylab="Cp",type='l')
plot(fwd.subset.summary$bic,xlab="Number of Variables (forward)",ylab="BIC",type='l')

par(mfrow=c(2,2))
plot(bk.subset.summary$rss,xlab="Number of Variables (backward)",ylab="RSS",type="l")
plot(bk.subset.summary$adjr2,xlab="Number of Variables (backward)",ylab="Adjusted RSq",type="l")
plot(bk.subset.summary$cp,xlab="Number of Variables (backward)",ylab="Cp",type='l')
plot(bk.subset.summary$bic,xlab="Number of Variables (backward)",ylab="BIC",type='l')

par(mfrow=c(2,2))
plot(ex.subset.summary$rss,xlab="Number of Variables (exhaustive)",ylab="RSS",type="l")
plot(ex.subset.summary$adjr2,xlab="Number of Variables (exhaustive)",ylab="Adjusted RSq",type="l")
plot(ex.subset.summary$cp,xlab="Number of Variables (exhaustive)",ylab="Cp",type='l')
plot(ex.subset.summary$bic,xlab="Number of Variables (exhaustive)",ylab="BIC",type='l')

Model Equations

Here we print the final equations for each model. Not, they are all the same.

for (mod.obj in model.list) {
  mod.name <- mod.obj[[1]]
  best.bic <- min(mod.obj[[3]]$bic)
  mod.num <- which.min(mod.obj[[3]]$bic)
  mod.cc <- coef(mod.obj[[2]], mod.num)
  
  mod.equation.format <- paste("Y =", paste(round(mod.cc[1],2), 
                      paste(round(mod.cc[-1],2), 
                      names(mod.cc[-1]), 
                      sep=" * ", collapse=" + "), 
                      sep=" + "), "+ e")
  
  print(paste("Model Selection Method: ",mod.name))
  print(paste("Min BIC:", best.bic))
  print(paste("Model Number: ", mod.num))
  print(paste("Model Equation: ", mod.equation.format))
  print("")
  
}
[1] "Model Selection Method:  Forward"
[1] "Min BIC: -373.710213587368"
[1] "Model Number:  3"
[1] "Model Equation:  Y = -0.27 + 1.27 * ShelveLocGood + 0.49 * CompPrice + -0.76 * Price + e"
[1] ""
[1] "Model Selection Method:  Backward"
[1] "Min BIC: -373.710213587368"
[1] "Model Number:  3"
[1] "Model Equation:  Y = -0.27 + 1.27 * ShelveLocGood + 0.49 * CompPrice + -0.76 * Price + e"
[1] ""
[1] "Model Selection Method:  Exhaustive"
[1] "Min BIC: -373.710213587368"
[1] "Model Number:  3"
[1] "Model Equation:  Y = -0.27 + 1.27 * ShelveLocGood + 0.49 * CompPrice + -0.76 * Price + e"
[1] ""
LS0tCnRpdGxlOiAiQ2Fyc2VhdHMgUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyBJbnRyb2R1Y3Rpb24KSGVyZSBpcyBteSBhbmFseXNpcyBvZiBgZGF0YS9DYXJzZWF0c19vcmcuY3N2YC4gCgojIyBDb25maWd1cmF0aW9uCgpgYGB7cn0KbGlicmFyeShsZWFwcykKbGlicmFyeShzdHJpbmdyKQpsaWJyYXJ5KGNhcmV0KQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkoRGF0YUV4cGxvcmVyKQpsaWJyYXJ5KGRwbHlyKQpsaWJyYXJ5KGdnRXh0cmEpCmxpYnJhcnkoUkNvbG9yQnJld2VyKQpsaWJyYXJ5KHBsb3RseSkKbGlicmFyeShjb3JycGxvdCkKbGlicmFyeShodG1sdG9vbHMpCmxpYnJhcnkoTUFTUykKYGBgCgojIyBTeXN0ZW0gSW5mb3JtYXRpb24KRHVlIHRvIHRoZSBsYXJnZSBudW1iZXIgb2YgbGlicmFyaWVzIGluIHVzZSBJIGhhdmUgcHJvdmlkZWQgc3lzdGVtIGluZm9ybWF0aW9uLgoKYGBge3J9CnNlc3Npb25JbmZvKCkKc2FwcGx5KGMoJ3JlcHInLCAnSVJkaXNwbGF5JywgJ0lSa2VybmVsJyksIGZ1bmN0aW9uKHApIHBhc3RlKHBhY2thZ2VWZXJzaW9uKHApKSkKYGBgCgojIERhdGEKSSBsb2FkIHRoZSBkYXRhIGludG8gYHJgLCBhbmQgZHJvcCB0aGUgIklEIiBjb2x1bW4uIAoKYGBge3J9CmNhcnNlYXRzIDwtIHJlYWQuY3N2KCJkYXRhL0NhcnNlYXRzX29yZy5jc3YiLCBoZWFkZXIgPSBULCBzdHJpbmdzQXNGYWN0b3JzID0gVCkKCmRyb3BzIDwtIGMoIlgiKQpjYXJzZWF0cyA8LSBjYXJzZWF0c1sgLCAhKG5hbWVzKGNhcnNlYXRzKSAlaW4lIGRyb3BzKV0KaGVhZChjYXJzZWF0cywgMTApCmBgYAoKSGVyZSBJIGNyZWF0ZSB0d28gbmV3IGRhdGEgZnJhbWVzIHRvIG1hbmFnZSBudW1lcmljIGFuZCBjYXRlZ29yaWNhbCBkYXRhLgoKYGBge3J9CiMgZ2V0IHZlY3RvcnMgb2YgY29udGludW91cyBhbmQgY2F0ZWdvcmljYWwgY29scwpudW1zIDwtIGRwbHlyOjpzZWxlY3RfaWYoY2Fyc2VhdHMsIGlzLm51bWVyaWMpCmNhdHMgPC0gZHBseXI6OnNlbGVjdF9pZihjYXJzZWF0cywgaXMuZmFjdG9yKQoKbnVtc1tzYW1wbGUobnJvdyhudW1zKSwgMTApLCBdCgpjYXRzW3NhbXBsZShucm93KGNhdHMpLCAxMCksIF0KYGBgCgpMZXQncyBnZXQgc29tZSBxdWljayBzdW1tYXJpZXMgb2YgZWFjaDoKCmBgYHtyfQpwcmludCgnTnVtZXJpYyBTdW1tYXJpZXMnKQpzdW1tYXJ5KG51bXMpCnByaW50KCdDYXRlZ29yaWNhbCBTdW1tYXJpZXMnKQpzdW1tYXJ5KGNhdHMpCmBgYAoKYGBge3J9CnN0cihudW1zKQpzdHIoY2F0cykKYGBgCgoKIyMgRGF0YSBEaW1lbnNpb25hbGl0eQpUaGlzIGNvbW1hbmQgaXMgdG8gaW5zcGVjdCB0aGUgZGlmZmVyZW50IGRhdGEgdHlwZXMgaW4gdGhlIGRhdGEuIAoKYGBge3J9CnN0cihjYXJzZWF0cykKcHJpbnQoIiIpCnByaW50KHBhc3RlKCdOdW1iZXIgb2YgQ29sdW1uczonLCBuY29sKGNhcnNlYXRzKSkpCnByaW50KHBhc3RlKCdOdW1iZXIgb2YgTnVtZXJpYyBDb2x1bW5zOicsIG5jb2wobnVtcykpKQpwcmludChwYXN0ZSgnTnVtYmVyIG9mIENhdGVnb3JpY2FsIENvbHVtbnM6JywgbmNvbChjYXRzKSkpCmBgYAoKYGBge3J9CmRpbShjYXJzZWF0cykKZGltKG51bXMpCmRpbShjYXRzKQpgYGAKCgpIZXJlJ3MgYSBxdWljayB3YXkgdG8gZXhhbWluZSBnZW5lcmFsIHByb3BlcnRpZXMgb2YgdGhlIGRhdGE6CgpgYGB7cn0KRGF0YUV4cGxvcmVyOjppbnRyb2R1Y2UoZGF0YT1jYXJzZWF0cykKYGBgCgpGaW5hbGx5LCBJIHdhbnQgdG8gbG9vayBhdCB0aGUgZmlyc3QgYW5kIGxhc3Qgcm93cyBvZiB0aGUgZGF0YSBzZXQuIEp1c3QgdG8gYmUgc2FmZToKCmBgYHtyfQpoZWFkKGNhcnNlYXRzLCAyKQp0YWlsKGNhcnNlYXRzLCAyKQoKaGVhZChudW1zLCAyKQp0YWlsKG51bXMsIDIpCgpoZWFkKGNhdHMsIDIpCnRhaWwoY2F0cywgMikKYGBgCgoKCiMgTnVtZXJpYyBQbG90dGluZwpJIHN0YXJ0IG91dCB3aXRoIGEgZmV3IGdlbmVyYWwgc2NhdHRlciBwbG90cy4KCmBgYHtyfQpwbG90X2x5KGRhdGE9Y2Fyc2VhdHMsIAogICAgICAgIHg9fkFnZSwgCiAgICAgICAgeT1+U2FsZXMsIAogICAgICAgIG1vZGUgPSAnbWFya2VycycsIAogICAgICAgIHR5cGUgPSAnc2NhdHRlcicsIAogICAgICAgIGNvbG9yPX5TaGVsdmVMb2MpICU+JQogICAgICAgICAgICBsYXlvdXQodGl0bGUgPSAiQWdlLCBTaGVsZiBMb2NhdGlvbiwgYW5kIFNhbGVzIFNjYXR0ZXIgUGxvdCIsIHdpZHRoPTkwMCkgCmBgYAoKVGhpcyBwbG90IGJlbG93IHNob3dzIGdvb2Qgc2VwYXJhdGlvbiBhbmQgYSB3ZWFrIGxpbmVhciB0cmVuZC4gVGhlc2UgdmFyaWFibGVzIGFyZSB3b3J0aCBpbnZlc3RpZ2F0aW5nLgoKYGBge3J9CnBsb3RfbHkoZGF0YT1jYXJzZWF0cywgCiAgICAgICAgeD1+UHJpY2UsIAogICAgICAgIHk9flNhbGVzLCAKICAgICAgICBtb2RlID0gJ21hcmtlcnMnLCAKICAgICAgICB0eXBlID0gJ3NjYXR0ZXInLCAKICAgICAgICBjb2xvcj1+U2hlbHZlTG9jKSAlPiUKICAgICAgICAgICAgbGF5b3V0KHRpdGxlID0gIlByaWNlLCBTaGVsZiBMb2NhdGlvbiwgYW5kIFNhbGVzIFNjYXR0ZXIgUGxvdCIsIHdpZHRoPTkwMCkKYGBgCgpIZXJlIHdlIGluc3BlY3QgdGhlIGRlbnNpdHkgb2YgdGhlICdQcmljZSB2cyBTYWxlcycgcmVsYXRpb25zaGlwOgoKYGBge3J9CnBsb3RfbHkoZGF0YT1jYXJzZWF0cywgCiAgICAgICAgeD1+UHJpY2UsIAogICAgICAgIHk9flNhbGVzLCAKICAgICAgICBtb2RlID0gJ21hcmtlcnMnLCAKICAgICAgICBzaXplID0gflByaWNlLAogICAgICAgIHR5cGUgPSAnc2NhdHRlcicsCiAgICAgICAgY29sb3JzID0gIkRhcmsyIiwKICAgICAgICBhbHBoYSA9IC42KSAlPiUKICAgICAgICAgICAgbGF5b3V0KHRpdGxlID0gIlByaWNlLCBVUywgYW5kIFNhbGVzIFNjYXR0ZXIgUGxvdCIsIHdpZHRoPTkwMCkgCmBgYAoKIyBOb3JtYWxpemF0aW9uCkkgY2hvb3NlIHRvIG5vcm1hbGl6ZSB0aGUgbnVtZXJpYyBkYXRhIGluIG9yZGVyIHRvIGJlIGFibGUgdG8gcGxvdCBlYWNoIHZhcmlhYmxlIG9uIHRoZSBzYW1lIHNjYWxlLiBUaGlzIHdpbGwgYWxsb3cgbWUgdG8gaW52ZXN0aWdhdGUgdGhlIHZhcmlhdGlvbiBvZiBlYWNoIHByZWRpY3RvciByZWxhdGl2ZSB0byBTYWxlcy4KCmBgYHtyfQpwcmVPYmogPC0gcHJlUHJvY2VzcyhudW1zLCBtZXRob2Q9YygiY2VudGVyIiwgInNjYWxlIikpCnNjYWxlZC5udW1zIDwtIHByZWRpY3QocHJlT2JqLCBudW1zKQoKaGVhZChzY2FsZWQubnVtcywgMikKdGFpbChzY2FsZWQubnVtcywgMikKYGBgCgpgYGB7cn0Kc3RyKHNjYWxlZC5udW1zKQpwcmludCgiIikKc3VtbWFyeShzY2FsZWQubnVtcykKYGBgCgojIyBEaXN0cmlidXRpb25zCkhlcmUgYXJlIHNjYWxlZCBkaXN0cmlidXRpb25zIChoaXN0b2dyYW1zIGFuZCBkZW5zaXR5IHBsb3RzKSBmb3IgZWFjaCBudW1lcmljIHZhcmlhYmxlLCBpbmNsdWRpbmcgU2FsZXMuIFRoZSB2YXJpYWJsZXMgcmVsYXRpbmcgdG8gbW9uZXkgKCQpIHRlbmQgdG8gYmUgYXBwcm94aW1hdGVseSBub3JtYWwuIE1hbnkgb3RoZXIgdmFyaWFibGVzIHRlbmQgdG8gYmUgYXBwcm94aW1hdGVseSB1bmlmb3JtLCB3aGljaCBkb2VzIG5vdCBib2RlIHdlbGwgZm9yIHRoZWlyIHByZWRpY3RpdmUgcG93ZXIuCgpgYGB7cn0Kc2NhbGVkLm51bXMgJT4lCiAgICB0aWR5cjo6Z2F0aGVyKCkgJT4lIAogICAgICAgIGdncGxvdChhZXMoeD12YWx1ZSx5PS4uZGVuc2l0eS4uKSkrCgogICAgICAgICAgICBnZ3RpdGxlKCdEaXN0cmlidXRpb25zIG9mIENvbnRpbm91cyBWYXJpYWJsZXMgKHNjYWxlZCknKSArCgogICAgICAgICAgICBmYWNldF93cmFwKH4ga2V5LCBzY2FsZXMgPSAiZnJlZSIpICsKICAgICAgICAgICAgZ2VvbV9oaXN0b2dyYW0oZmlsbD1JKCJvcmFuZ2UiKSwgY29sPUkoImJsYWNrIiksIGJpbnMgPSA1MCkgKwoKICAgICAgICAgICAgZmFjZXRfd3JhcCh+IGtleSwgc2NhbGVzID0gImZyZWUiKSArCiAgICAgICAgICAgIGdlb21fZGVuc2l0eShjb2xvcj0iYmx1ZSIsIGZpbGw9J2xpZ2h0IGJsdWUnLCBhbHBoYSA9IDAuNCkKYGBgCgpIZXJlIHdlIHBsb3QgYWxsIG51bWVyaWMgdmFyaWFibGVzIGFnYWluc3QgdGhlaXIgZGlzdHJpYnV0aW9ucy4gVGhpcyBpcyBqdXN0IGFub3RoZXIgd2F5IHRvIGV4YW1pbmUgdGhlIGluZm9ybWF0aW9uIHNob3duIGFib3ZlLgoKYGBge3J9CnNjYWxlZC5udW1zICU+JQogICAgdGlkeXI6OmdhdGhlcigpICU+JSAKICAgICAgICBwbG90X2x5KHg9fmtleSwgeT1+dmFsdWUsCiAgICAgICAgICAgICAgICB0eXBlID0gImJveCIsIAogICAgICAgICAgICAgICAgYm94cG9pbnRzID0gImFsbCIsIAogICAgICAgICAgICAgICAgaml0dGVyID0gMC40LAogICAgICAgICAgICAgICAgcG9pbnRwb3MgPSAwLAogICAgICAgICAgICAgICAgY29sb3IgPSB+a2V5LAogICAgICAgICAgICAgICAgY29sb3JzID0gIkRhcmsyIikgJT4lCiAgICAgICAgICAgICAgICAgICAgICBzdWJwbG90KHNoYXJlWCA9IFRSVUUpICAlPiUKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGxheW91dCh0aXRsZSA9ICJOdW1lcmljIFZhcmlhYmxlIERpc3RyaWJ1dGlvbnMgKHNjYWxlZCkiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzPWxpc3QodGl0bGU9J1N0YW5kYXJkIERldmlhdGlvbicpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeGF4aXM9bGlzdCh0aXRsZT0nVmFyaWFibGUnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGF1dG9zaXplPUZBTFNFLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgd2lkdGg9OTAwLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgaGVpZ2h0PTUwMCkKYGBgCgojIyBTY2F0dGVycGxvdHMKCkhlcmUgd2UgcGxvdCBhbGwgbnVtZXJpYyB2YXJpYWJsZXMgYWdhaW5zdCBTYWxlcyAoc2NhbGVkKS4gVGhpcyBhbGxvd3MgdXMgdG8gaW52ZXN0aWdhdGUgcG9zc2libGUgbGluZWFyIHJlbGF0aW9uc2hpcHMgYmV0d2VlbiB0aGF0IHZhcmlhYmxlIGFuZCBTYWxlcy4gQXMgc2hvd24gYmVsb3csIG9ubHkgJ1ByaWNlJyBhcHBlYXJzIHRvIGhhdmUgYSBsaW5lYXIgcmVsYXRpb25zaGlwIHdvcnRoIGludmVzdGlnYXRpbmcuIFRoaXMgdG9vayBtZSBzbyBsb25nIHRvIGZpZ3VyZSBvdXQuIAoKYGBge3J9CgpudW1lcmljLnNjYXR0ZXJwbG90cyA8LSBodG1sdG9vbHM6OnRhZ0xpc3QoKQpjb3VudCA9IDEKCmZvciAoaSBpbiBuYW1lcyhzY2FsZWQubnVtc1ssLTFdKSkgewogIAogIG51bWVyaWMuc2NhdHRlcnBsb3RzW1tjb3VudF1dIDwtIHBsb3RfbHkoc2NhbGVkLm51bXMsIHg9c2NhbGVkLm51bXNbLGldLCB5PXNjYWxlZC5udW1zJFNhbGVzLAogICAgICAgICAgICAgICAgICAgIGNvbG9ycyA9ICdSZFlsR24nLAogICAgICAgICAgICAgICAgICAgIG1vZGUgPSAnbWFya2VycycsCiAgICAgICAgICAgICAgICAgICAgdHlwZSA9ICdzY2F0dGVyJywKICAgICAgICAgICAgICAgICAgICBzaXplID0gc2NhbGVkLm51bXMkU2FsZXNeMiwKICAgICAgICAgICAgICAgICAgICBjb2xvciA9IHNjYWxlZC5udW1zJFNhbGVzLAogICAgICAgICAgICAgICAgICAgIG1hcmtlciA9IGxpc3QobGluZSA9IGxpc3QoY29sb3IgPSAnYmxhY2snLHdpZHRoID0gMikpLAogICAgICAgICAgICAgICAgICAgIG5hbWU9cGFzdGUoaSkpICU+JQogICAgICAgICAgbGF5b3V0KHRpdGxlID0gcGFzdGUoaSwgInZzIFNhbGVzIChzY2FsZWQpPGJyPlNpemU9U2FsZXNeMiIpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeWF4aXM9bGlzdCh0aXRsZT0nU2FsZXMnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHhheGlzPWxpc3QodGl0bGU9aSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzaG93bGVnZW5kID0gRkFMU0UpCiAgCiAgY291bnQgPSBjb3VudCArIDEKICAKfQoKbnVtZXJpYy5zY2F0dGVycGxvdHMKCmBgYAoKQnkgYWRkaW5nIG5haXZlIHJlZ3Jlc3Npb24gbGluZXMgdG8gYSBmZXcgc2NhdHRlciBwbG90cyB3ZSBjYW4gY29uZmlybSBvdXIgc3VzcGljaW9uczoKCmBgYHtyfQpmaXQuUG9wIDwtIGxtKFNhbGVzIH4gUG9wdWxhdGlvbiwgZGF0YSA9IHNjYWxlZC5udW1zKQpmaXQuQWdlIDwtIGxtKFNhbGVzIH4gQWdlLCBkYXRhID0gc2NhbGVkLm51bXMpCmZpdC5Db21wUHJpY2UgPC0gbG0oU2FsZXMgfiBDb21wUHJpY2UsIGRhdGEgPSBzY2FsZWQubnVtcykKZml0LlByaWNlIDwtIGxtKFNhbGVzIH4gUHJpY2UsIGRhdGEgPSBzY2FsZWQubnVtcykKCnJlZ3Jlc3Npb24uc2NhdHRlcnBsb3RzIDwtIGh0bWx0b29sczo6dGFnTGlzdCgpCgpyZWdyZXNzaW9uLnNjYXR0ZXJwbG90c1tbMV1dIDwtICBwbG90X2x5KHNjYWxlZC5udW1zLCAKICAgICAgICAgIHggPSB+UG9wdWxhdGlvbiwKICAgICAgICAgIG5hbWUgPSAnUG9wdWxhdGlvbiB2cyBTYWxlcyBSZWdyZXNzaW9uIExpbmUnKSAlPiUgCiAgICAgICAgICAgICAgYWRkX21hcmtlcnMoeSA9IH5TYWxlcywKICAgICAgICAgICAgICAgICAgICBuYW1lID0gJ1BvcHVsYXRpb24gdnMgU2FsZXMgT2JzZXJ2YXRpb25zJykgJT4lIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgYWRkX2xpbmVzKHggPSB+UG9wdWxhdGlvbiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5ID0gZml0dGVkKGZpdC5Qb3ApKSAlPiUKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGxheW91dCh0aXRsZSA9ICJQb3B1bGF0aW9uIHZzIFNhbGVzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5YXhpcz1saXN0KHRpdGxlPSdTYWxlcycsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeGF4aXM9bGlzdCh0aXRsZT0nUG9wdWxhdGlvbicpKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzaG93bGVnZW5kID0gRkFMU0UpCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgICAgICAgICAKCnJlZ3Jlc3Npb24uc2NhdHRlcnBsb3RzW1syXV0gPC0gIHBsb3RfbHkoc2NhbGVkLm51bXMsIAogICAgICAgICAgeCA9IH5BZ2UsCiAgICAgICAgICBuYW1lID0gJ0FnZSB2cyBTYWxlcyBSZWdyZXNzaW9uIExpbmUnKSAlPiUgCiAgICAgICAgICAgICAgYWRkX21hcmtlcnMoeSA9IH5TYWxlcywKICAgICAgICAgICAgICAgICAgICBuYW1lID0gJ0FnZSB2cyBTYWxlcyBPYnNlcnZhdGlvbnMnKSAlPiUgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBhZGRfbGluZXMoeCA9IH5BZ2UsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeSA9IGZpdHRlZChmaXQuQWdlKSkgJT4lCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBsYXlvdXQodGl0bGUgPSAiQWdlIHZzIFNhbGVzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5YXhpcz1saXN0KHRpdGxlPSdTYWxlcycsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeGF4aXM9bGlzdCh0aXRsZT0nQWdlJykpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHNob3dsZWdlbmQgPSBGQUxTRSkKCnJlZ3Jlc3Npb24uc2NhdHRlcnBsb3RzW1szXV0gPC0gIHBsb3RfbHkoc2NhbGVkLm51bXMsIAogICAgICAgICAgeCA9IH5Db21wUHJpY2UsCiAgICAgICAgICBuYW1lID0gJ0NvbXBQcmljZSB2cyBTYWxlcyBSZWdyZXNzaW9uIExpbmUnKSAlPiUgCiAgICAgICAgICAgICAgYWRkX21hcmtlcnMoeSA9IH5TYWxlcywKICAgICAgICAgICAgICAgICAgICBuYW1lID0gJ0NvbXBQcmljZSB2cyBTYWxlcyBPYnNlcnZhdGlvbnMnKSAlPiUgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBhZGRfbGluZXMoeCA9IH5Db21wUHJpY2UsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeSA9IGZpdHRlZChmaXQuQ29tcFByaWNlKSkgJT4lCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBsYXlvdXQodGl0bGUgPSAiQ29tcFByaWNlIHZzIFNhbGVzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5YXhpcz1saXN0KHRpdGxlPSdTYWxlcycsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeGF4aXM9bGlzdCh0aXRsZT0nQ29tcFByaWNlJykpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHNob3dsZWdlbmQgPSBGQUxTRSkKCnJlZ3Jlc3Npb24uc2NhdHRlcnBsb3RzW1s0XV0gPC0gIHBsb3RfbHkoc2NhbGVkLm51bXMsIAogICAgICAgICAgeCA9IH5QcmljZSwKICAgICAgICAgIG5hbWUgPSAnUHJpY2UgdnMgU2FsZXMgUmVncmVzc2lvbiBMaW5lJykgJT4lIAogICAgICAgICAgICAgIGFkZF9tYXJrZXJzKHkgPSB+U2FsZXMsCiAgICAgICAgICAgICAgICAgICAgbmFtZSA9ICdQcmljZSB2cyBTYWxlcyBPYnNlcnZhdGlvbnMnKSAlPiUgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBhZGRfbGluZXMoeCA9IH5QcmljZSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5ID0gZml0dGVkKGZpdC5QcmljZSkpICU+JQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbGF5b3V0KHRpdGxlID0gIlByaWNlIHZzIFNhbGVzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5YXhpcz1saXN0KHRpdGxlPSdTYWxlcycsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeGF4aXM9bGlzdCh0aXRsZT0nUHJpY2UnKSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc2hvd2xlZ2VuZCA9IEZBTFNFKQoKCnJlZ3Jlc3Npb24uc2NhdHRlcnBsb3RzCmBgYAoKTGV0J3MgY29tcGFyZSB0aGUgc2xvcGVzOgoKYGBge3J9CnNjYWxlZC5udW1zICU+JQogICAgcGxvdF9seSh5ID0gflNhbGVzKSAlPiUKICAgICAgYWRkX2xpbmVzKHg9IH5Qb3B1bGF0aW9uLCB5ID0gZml0dGVkKGZpdC5Qb3ApLCAKICAgICAgICAgICAgICAgIG5hbWUgPSAiZml0LlBvcCBzbG9wZSIsIGxpbmUgPSBsaXN0KHNoYXBlID0gImxpbmVhciIpKSAlPiUKICAgICAgYWRkX2xpbmVzKHg9IH5BZ2UsIHkgPSBmaXR0ZWQoZml0LkFnZSksIAogICAgICAgICAgICAgICAgbmFtZSA9ICJmaXQuQWdlIHNsb3BlIiwgbGluZSA9IGxpc3Qoc2hhcGUgPSAibGluZWFyIikpICU+JQogICAgICBhZGRfbGluZXMoeD0gfkNvbXBQcmljZSwgeSA9IGZpdHRlZChmaXQuQ29tcFByaWNlKSwgCiAgICAgICAgICAgICAgICBuYW1lID0gImZpdC5Db21wUHJpY2Ugc2xvcGUiLCBsaW5lID0gbGlzdChzaGFwZSA9ICJsaW5lYXIiKSkgJT4lCiAgICAgIGFkZF9saW5lcyh4PSB+UHJpY2UsIHkgPSBmaXR0ZWQoZml0LlByaWNlKSwgCiAgICAgICAgICAgICAgICBuYW1lID0gImZpdC5QcmljZSBzbG9wZSIsIGxpbmUgPSBsaXN0KHNoYXBlID0gImxpbmVhciIpKSAlPiUKICAgICAgICAgICAgCiAgICBsYXlvdXQodGl0bGUgPSAiUmVncmVzc2lvbiBMaW5lcyB2cyBTYWxlcyAoc2NhbGVkKSIsCiAgICAgICAgICAgYXV0b3NpemU9RkFMU0UsCiAgICAgICAgICAgd2lkdGg9OTAwLAogICAgICAgICAgIHlheGlzPWxpc3QodGl0bGU9J1NhbGVzJyksCiAgICAgICAgICAgeGF4aXM9bGlzdCh0aXRsZT0nU2NhbGVkIE51bWVyaWMgVmFyaWFibGUnKSkKYGBgCgpIZXJlJ3MgYSBwcmV0dHkgZ3JhcGhpYyB0aGF0IGRvZXNuJ3QgaGVscCBtZSB1bmRlcnN0YW5kIGFueXRoaW5nIGFib3V0IHRoZSBkYXRhLgoKYGBge3J9CnkgPSBzY2FsZWQubnVtcyRTYWxlcwp4ID0gc2NhbGVkLm51bXMkUHJpY2UKCnMgPC0gc3VicGxvdCgKICBwbG90X2x5KHggPSB4LCBjb2xvciA9IEkoImJsYWNrIiksIHR5cGUgPSAnaGlzdG9ncmFtJyksIAogIHBsb3RseV9lbXB0eSgpLCAKICBwbG90X2x5KHggPSB4LCB5ID0geSwgdHlwZSA9ICdoaXN0b2dyYW0yZGNvbnRvdXInLCBzaG93c2NhbGUgPSBGKSwgCiAgcGxvdF9seSh5ID0geSwgY29sb3IgPSBJKCJibGFjayIpLCB0eXBlID0gJ2hpc3RvZ3JhbScpLAogIG5yb3dzID0gMiwgaGVpZ2h0cyA9IGMoMC4yLCAwLjgpLCB3aWR0aHMgPSBjKDAuOCwgMC4yKSwKICBzaGFyZVggPSBUUlVFLCBzaGFyZVkgPSBUUlVFLCB0aXRsZVggPSBGQUxTRSwgdGl0bGVZID0gRkFMU0UpCgpsYXlvdXQocywgc2hvd2xlZ2VuZCA9IEZBTFNFLCBhdXRvc2l6ZT1GQUxTRSwKICAgICAgICAgICB3aWR0aD03MDAsCiAgICAgICAgICAgaGVpZ2h0PTUwMCwgCiAgICAgICAgICAgeWF4aXM9bGlzdCh0aXRsZT0nU2FsZXMnKSwKICAgICAgICAgICB4YXhpcz1saXN0KHRpdGxlPSdQcmljZScpKQpgYGAKCiMgQ2F0ZWdvcmljYWwgUGxvdHRpbmcKRmlyc3QsIGxldCdzIGNyZWF0ZSBhIGRhdGEgZnJhbWUgdGhhdCB3ZSBjYW4gdXNlOgpgYGB7cn0KY2F0ZWdvcmljYWwuYnkuc2FsZXMgPSBjYmluZChTYWxlcyA9IHNjYWxlZC5udW1zJFNhbGVzLCBjYXRzKQpzdHIoY2F0ZWdvcmljYWwuYnkuc2FsZXMpCmBgYAoKSGVyZSB3ZSBjYW4gc2VlIGFsbCBjYXRlZ29yaWNhbCBieSBTYWxlcy4gV2Ugc3VzcGVjdGVkIHRoYXQgJ1NoZWx2ZUxvYycgd291bGQgYmUgaW1wb3J0YW50IGJhc2VkIG9uIG9uZSBvZiB0aGUgZWFybHkgc2NhdHRlciBwbG90cy4gSXQgc2VlbXMgdGhhdCB0aGlzIGlzIHRoZSBjYXNlLgoKYGBge3J9CgpjYXRlZ29yaWNhbC5ib3hwbG90cyA8LSBodG1sdG9vbHM6OnRhZ0xpc3QoKQpjb3VudCA9IDEKCmZvciAoaSBpbiBuYW1lcyhjYXRlZ29yaWNhbC5ieS5zYWxlc1ssLTFdKSkgewogIAogIGNhdGVnb3JpY2FsLmJveHBsb3RzW1tjb3VudF1dIDwtIHBsb3RfbHkoY2F0ZWdvcmljYWwuYnkuc2FsZXMsIHg9Y2F0ZWdvcmljYWwuYnkuc2FsZXNbLGldLCB5PWNhdGVnb3JpY2FsLmJ5LnNhbGVzJFNhbGVzLAogICAgICAgICAgICAgICAgICAgICAgICB0eXBlID0gImJveCIsIAogICAgICAgICAgICAgICAgICAgICAgICBib3hwb2ludHMgPSAiYWxsIiwKICAgICAgICAgICAgICAgICAgICAgICAgaml0dGVyID0gLjIsCiAgICAgICAgICAgICAgICAgICAgICAgIHBvaW50cG9zID0gMCwKICAgICAgICAgICAgICAgICAgICAgICAgY29sb3IgPWNhdGVnb3JpY2FsLmJ5LnNhbGVzWyxpXSwKICAgICAgICAgICAgICAgICAgICAgICAgY29sb3JzPSdTZXQxJywgCiAgICAgICAgICAgICAgICAgICAgICAgIG5hbWU9cGFzdGUoaSkpICU+JQogICAgICAgICAgICAgICAgICAgICAgICAgICAgbGF5b3V0KHRpdGxlID0gcGFzdGUoaSwgInZzIFNhbGVzIChzY2FsZWQpIiksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzaG93bGVnZW5kID0gVFJVRSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzPWxpc3QodGl0bGU9J1NhbGVzIFN0YW5kYXJkIERldmlhdGlvbicpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeGF4aXM9bGlzdCh0aXRsZT1pKSkKICAKICBjb3VudCA9IGNvdW50ICsgMQogIAp9CgpjYXRlZ29yaWNhbC5ib3hwbG90cwpgYGAKCkhlcmUncyB0aGUgc2FtZSB0aGluZywgYnV0IG1vcmUgbXVzaWNhbGx5OgoKYGBge3J9CgpjYXRlZ29yaWNhbC52aW9saW5zIDwtIGh0bWx0b29sczo6dGFnTGlzdCgpCmNvdW50ID0gMQoKZm9yIChpIGluIG5hbWVzKGNhdGVnb3JpY2FsLmJ5LnNhbGVzWywtMV0pKSB7CiAgCiAgY2F0ZWdvcmljYWwudmlvbGluc1tbY291bnRdXSA8LSBwbG90X2x5KGNhdGVnb3JpY2FsLmJ5LnNhbGVzLCB4PWNhdGVnb3JpY2FsLmJ5LnNhbGVzWyxpXSwgeT1jYXRlZ29yaWNhbC5ieS5zYWxlcyRTYWxlcywKICAgICAgICAgICAgICAgICAgICAgICAgc3BsaXQgPSBjYXRlZ29yaWNhbC5ieS5zYWxlc1ssaV0sCiAgICAgICAgICAgICAgICAgICAgICAgIHR5cGUgPSAndmlvbGluJywKICAgICAgICAgICAgICAgICAgICAgICAgY29sb3JzPSdTZXQxJywgCiAgICAgICAgICAgICAgICAgICAgICAgIG5hbWU9cGFzdGUoaSksCiAgICAgICAgICAgICAgICAgICAgICAgIGJveCA9IGxpc3QodmlzaWJsZSA9IFRSVUUpLAogICAgICAgICAgICAgICAgICAgICAgICBtZWFubGluZSA9IGxpc3QodmlzaWJsZSA9IFRSVUUpKSAlPiUgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBsYXlvdXQoeGF4aXMgPSBsaXN0KHRpdGxlID0gIlVTIiksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeWF4aXMgPSBsaXN0KHRpdGxlID0gIlNhbGVzIix6ZXJvbGluZSA9IEZBTFNFKSkKCiAgCiAgY291bnQgPSBjb3VudCArIDEKICAKfQoKY2F0ZWdvcmljYWwudmlvbGlucwpgYGAKCiMgTGluZWFyIFJlZ3Jlc3Npb24KCkZpcnN0LCBsZXQncyBtZXJnZSB0aGUgZGF0YSBzZXQgaW50byBhIHNpbmdsZSBkYXRhIGZyYW1lCmBgYHtyfQoKc2NhbGVkLm1lcmdlZCA8LSBjYmluZChjYXRlZ29yaWNhbC5ieS5zYWxlc1ssLTFdLCBzY2FsZWQubnVtcykKc3RyKHNjYWxlZC5tZXJnZWQpCmBgYAoKYGBge3J9CmhlYWQobnVtcywgMikKdGFpbChudW1zLCAyKQoKaGVhZChzY2FsZWQubWVyZ2VkLCAyKQp0YWlsKHNjYWxlZC5tZXJnZWQsIDIpCmBgYAoKCkZpcnN0LCBsZXQncyBsb29rIGF0IHNvbWUgdGhpbmdzIHRoYXQgbWF5IGdpdmUgdXMgdHJvdWJsZS4gTHVja2lseSBpdCBsb29rcyBsaWtlIHRoZSBvbmx5IHNlcmlvdXMgY29ycmVsYXRpb24gaXMgd2l0aCBvdXIgZGVwZW5kZW50IHZhcmlhYmxlLiBXZSdsbCB3YW50IHRvIHdhdGNoIHRoZSAnUHJpY2UnIHZzICdDb21wUHJpY2UnIHJlbGF0aW9uc2hpcC4KCmBgYHtyfQpyZXMgPC0gY29yKHNjYWxlZC5udW1zKQoKY29ycnBsb3QocmVzLCB0eXBlID0gInVwcGVyIiwgb3JkZXIgPSAiaGNsdXN0IiwgCiAgICAgICAgIHRsLmNvbCA9ICJibGFjayIsIHRsLnNydCA9IDQ1KQpgYGAKCkl0IGFwcGVhcnMgdGhhdCByZXNpZHVhbHMgYXJlIHJvdWdobHkgc3ltbWV0cmljYWwgYXJvdW5kIDAuIFRoYXQncyBzdHJhbmdlLiBNb3N0bHkgZHVlIHRvIGEgcmVsYXRpdmVseSBwb29yIG92ZXJhbGwgZml0LiBOb3RlIGhvdyBjbG9zZSB0byB6ZXJvIG1vc3Qgb2YgdGhlIGNvZWZmaWNpZW50IGVzdGltYXRlcyBhcmUuCgpgYGB7cn0Kc2ltcGxlLmxtIDwtIGxtKFNhbGVzfi4sIGRhdGE9c2NhbGVkLm1lcmdlZCkKc2ltcGxlLnN1bW1hcnkgPC0gc3VtbWFyeShzaW1wbGUubG0pCgpwcmludChzaW1wbGUuc3VtbWFyeSkKYGBgCgojIyBMaW5lYXIgTW9kZWxzIGFuZCBTdWJzZXRzCgpMZXQncyBkbyB0aGUgc2FtZSB0aGluZywgYnV0IGNvbnRyb2wgdGhlIHN1YnNldHMgdXNpbmcgYGxlYXBzYAoKYGBge3J9CnJlZ2ZpdC5mdWxsPXJlZ3N1YnNldHMoU2FsZXN+LiwgZGF0YT1zY2FsZWQubWVyZ2VkLCBudm1heD01KQpyZWcuc3VtbWFyeT1zdW1tYXJ5KHJlZ2ZpdC5mdWxsKQoKcHJpbnQocmVnLnN1bW1hcnkpCmBgYAoKV2UnbGwganVzdCB0YWtlIGNvZGUgc3RyYWlnaHQgZnJvbSB0aGUgZXhhbXBsZSBvbiBDYW52YXMuLi4KCmBgYHtyfQpwYXIobWZyb3c9YygyLDIpKQpwbG90KHJlZy5zdW1tYXJ5JHJzcyx4bGFiPSJOdW1iZXIgb2YgVmFyaWFibGVzIix5bGFiPSJSU1MiLHR5cGU9ImwiKQpwbG90KHJlZy5zdW1tYXJ5JGFkanIyLHhsYWI9Ik51bWJlciBvZiBWYXJpYWJsZXMiLHlsYWI9IkFkanVzdGVkIFJTcSIsdHlwZT0ibCIpCnBsb3QocmVnLnN1bW1hcnkkY3AseGxhYj0iTnVtYmVyIG9mIFZhcmlhYmxlcyIseWxhYj0iQ3AiLHR5cGU9J2wnKQpwbG90KHJlZy5zdW1tYXJ5JGJpYyx4bGFiPSJOdW1iZXIgb2YgVmFyaWFibGVzIix5bGFiPSJCSUMiLHR5cGU9J2wnKQpgYGAKCmBgYHtyfQpwbG90KHJlZ2ZpdC5mdWxsLHNjYWxlPSJyMiIpCnBsb3QocmVnZml0LmZ1bGwsc2NhbGU9ImFkanIyIikKcGxvdChyZWdmaXQuZnVsbCxzY2FsZT0iQ3AiKQpwbG90KHJlZ2ZpdC5mdWxsLHNjYWxlPSJiaWMiKQpgYGAKCgojIEludGVyYWN0aW9uIFRlcm1zCgpIZXJlIHdlIGRlZmluZSBhIG5ldyBtb2RlbCB3aXRoIHNvbWUgaW50ZXJhY3Rpb24gdGVybXM6IAogICAgYS4gIEluY29tZSBhbmQgQWR2ZXJ0aXNpbmcKICAgIGIuICBJbmNvbWUgYW5kIENvbXBQcmljZQogICAgYy4gICBQcmljZSBhbmQgQWdlCgoKYGBge3J9CmludGVyYWN0aW9uLmxtIDwtIGxtKFNhbGVzfi4gKyBJbmNvbWUqQWR2ZXJ0aXNpbmcgKyBJbmNvbWUqQ29tcFByaWNlICsgUHJpY2UqQWdlLCBkYXRhPXNjYWxlZC5tZXJnZWQpCmludGVyYWN0aW9uLnN1bW1hcnkgPC0gc3VtbWFyeShpbnRlcmFjdGlvbi5sbSkKCnByaW50KGludGVyYWN0aW9uLnN1bW1hcnkpCmBgYAoKYGBge3J9CmludGVyYWN0aW9uLmxtLnN1YnNldHMgPC0gcmVnc3Vic2V0cyhTYWxlc34uICsgSW5jb21lKkFkdmVydGlzaW5nICsgSW5jb21lKkNvbXBQcmljZSArIFByaWNlKkFnZSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBkYXRhPXNjYWxlZC5tZXJnZWQsIG52bWF4PTUpCgppbnRlcmFjdGlvbi5zdWJzZXRzLnN1bW1hcnkgPC0gc3VtbWFyeShpbnRlcmFjdGlvbi5sbS5zdWJzZXRzKQoKcHJpbnQoaW50ZXJhY3Rpb24uc3Vic2V0cy5zdW1tYXJ5KQpgYGAKCgojIyBWYXJpYWJsZSBTaWduaWZpY2FuY2UKQmVsb3cgd2UgcHJpbnQgdGhlIGNvZWZmaWNpZW50cyBmb3IgdGhlIDV0aCBtb2RlbCB1c2luZyB0aGUgZGVmYXVsdCBtb2RlbCBzZWxlY3Rpb24gY3JpdGVyaWEuIEFsbCBjb2VmZmljaWVudHMgYXJlIHJlbGF0aXZlbHkgc21hbGwsIGFzIHdlIHdvdWxkIGV4cGVjdCBmcm9tIHRoZSBFREEgYWJvdmUuIFRoaXMgcHJldHR5IG11Y2ggY29uZmlybXMgd2hhdCBJIHdvdWxkIGhhdmUgZ3Vlc3NlZCBieSBsb29raW5nIGF0IHRoZSBkYXRhIGFnYWluc3Qgc2FsZXMuIFdlIHN0aWxsIHdhbnQgdG8gd2F0Y2ggb3V0IGZvciBjb25mb3VuZGluZyBiZXR3ZWVuICdQcmljZScgYW5kICdDb21wUHJpY2UuJwoKCmBgYHtyfQpjb2VmKGludGVyYWN0aW9uLmxtLnN1YnNldHMsIDUpCmBgYAoKIyMgU2Vjb25kIEludGVyYWN0aW9uIE1vZGVsCgpGaXJzdCwgZHJvcCBjb2x1bW5zIHVubmVlZGVkIGZyb20gYW5hbHlzaXM6CgpgYGB7cn0Kc2NhbGVkLm1lcmdlZC5zbGltIDwtIHNjYWxlZC5tZXJnZWRbICwgLXdoaWNoKG5hbWVzKHNjYWxlZC5tZXJnZWQpICVpbiUgYygiVVMiLCJVcmJhbiIpKV0KYGBgCgoKQSBmZXcgaHlwZXIgcGFyYW1ldGVycyB3ZSdkIGxpa2UgdG8gYmUgY29uc2lzdGVudCBmb3IgYWxsIG1vZGVscwoKYGBge3J9Cm52bWF4IDwtIDMKYGBgCgojIyBGb3J3YXJkIFNlbGVjdGlvbjoKCmBgYHtyfQppbnRlcmFjdGlvbi5zdWJzZXQuZndkIDwtIHJlZ3N1YnNldHMoU2FsZXN+LiArIEluY29tZSpBZHZlcnRpc2luZyArIEluY29tZSpDb21wUHJpY2UgKyBJbmNvbWUqQWdlLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZGF0YT1zY2FsZWQubWVyZ2VkLnNsaW0sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbnZtYXggPSBudm1heCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1ldGhvZD0iZm9yd2FyZCIpCgpmd2Quc3Vic2V0LnN1bW1hcnkgPC0gc3VtbWFyeShpbnRlcmFjdGlvbi5zdWJzZXQuZndkKQoKY29lZihpbnRlcmFjdGlvbi5zdWJzZXQuZndkLCAxOm52bWF4KQpgYGAKCiMjIEJhY2t3YXJkIFNlbGVjdGlvbjoKClRoaXMgaXMgcmVhbGx5IHN0cmFuZ2UuIEkgY2FuJ3Qgc2VlbSB0byBmaW5kIGFueSBkb2N1bWVudGF0aW9uIGFib3V0IHRoaXMsIGJ1dCBpdCBhcHBlYXJzIHRoYXQgdGhpcyBtb2RlbCBpcyBhY3R1YWxseSAnZm9yd2FyZC4nCgpgYGB7cn0KaW50ZXJhY3Rpb24uc3Vic2V0LmJrIDwtIHJlZ3N1YnNldHMoU2FsZXN+LiArIEluY29tZSpBZHZlcnRpc2luZyArIEluY29tZSpDb21wUHJpY2UgKyBJbmNvbWUqQWdlLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZGF0YT1zY2FsZWQubWVyZ2VkLnNsaW0sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbnZtYXggPSBudm1heCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1ldGhvZD0iYmFja3dhcmQiKQoKYmsuc3Vic2V0LnN1bW1hcnkgPC0gc3VtbWFyeShpbnRlcmFjdGlvbi5zdWJzZXQuYmspCgpjb2VmKGludGVyYWN0aW9uLnN1YnNldC5iaywgMTpudm1heCkKYGBgCgojIyBFeGhhc3V0aXZlCgpgYGB7cn0KaW50ZXJhY3Rpb24uc3Vic2V0LmV4IDwtIHJlZ3N1YnNldHMoU2FsZXN+LiArIEluY29tZSpBZHZlcnRpc2luZyArIEluY29tZSpDb21wUHJpY2UgKyBJbmNvbWUqQWdlLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZGF0YT1zY2FsZWQubWVyZ2VkLnNsaW0sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbnZtYXggPSBudm1heCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1ldGhvZD0iZXhoYXVzdGl2ZSIpCgpleC5zdWJzZXQuc3VtbWFyeSA8LSBzdW1tYXJ5KGludGVyYWN0aW9uLnN1YnNldC5leCkKY29lZihpbnRlcmFjdGlvbi5zdWJzZXQuZXgsIDE6bnZtYXgpCmBgYAoKVGhpcyBpcyBhIGxpc3QgdGhhdCBtYXkgY29tZSBpbiBoYW5keS4KCmBgYHtyfQptb2RlbC5saXN0ID0gbGlzdChsaXN0KCJGb3J3YXJkIiwgaW50ZXJhY3Rpb24uc3Vic2V0LmZ3ZCwgZndkLnN1YnNldC5zdW1tYXJ5KSwgCiAgICAgICAgICAgICAgICAgIGxpc3QoIkJhY2t3YXJkIiwgaW50ZXJhY3Rpb24uc3Vic2V0LmJrLCBiay5zdWJzZXQuc3VtbWFyeSksIAogICAgICAgICAgICAgICAgICBsaXN0KCJFeGhhdXN0aXZlIiwgaW50ZXJhY3Rpb24uc3Vic2V0LmV4LCBleC5zdWJzZXQuc3VtbWFyeSkpCmBgYAoKIyBFdmFsdWF0aW9uIE1ldHJpYyBQbG90dGluZwoKSXQgaXMgaW50ZXJlc3RpbmcgdG8gc2VlIHRoYXQgZWFjaCBtb2RlbCBzZWxlY3RlZCB0aGUgc2FtZSB2YXJpYWJsZXMsIGluIHRoZSBzYW1lIG9yZGVyLgoKYGBge3J9CnBhcihtZnJvdz1jKDIsMikpCnBsb3QoZndkLnN1YnNldC5zdW1tYXJ5JHJzcyx4bGFiPSJOdW1iZXIgb2YgVmFyaWFibGVzIChmb3J3YXJkKSIseWxhYj0iUlNTIix0eXBlPSJsIikKcGxvdChmd2Quc3Vic2V0LnN1bW1hcnkkYWRqcjIseGxhYj0iTnVtYmVyIG9mIFZhcmlhYmxlcyAoZm9yd2FyZCkiLHlsYWI9IkFkanVzdGVkIFJTcSIsdHlwZT0ibCIpCnBsb3QoZndkLnN1YnNldC5zdW1tYXJ5JGNwLHhsYWI9Ik51bWJlciBvZiBWYXJpYWJsZXMgKGZvcndhcmQpIix5bGFiPSJDcCIsdHlwZT0nbCcpCnBsb3QoZndkLnN1YnNldC5zdW1tYXJ5JGJpYyx4bGFiPSJOdW1iZXIgb2YgVmFyaWFibGVzIChmb3J3YXJkKSIseWxhYj0iQklDIix0eXBlPSdsJykKYGBgCgpgYGB7cn0KcGFyKG1mcm93PWMoMiwyKSkKcGxvdChiay5zdWJzZXQuc3VtbWFyeSRyc3MseGxhYj0iTnVtYmVyIG9mIFZhcmlhYmxlcyAoYmFja3dhcmQpIix5bGFiPSJSU1MiLHR5cGU9ImwiKQpwbG90KGJrLnN1YnNldC5zdW1tYXJ5JGFkanIyLHhsYWI9Ik51bWJlciBvZiBWYXJpYWJsZXMgKGJhY2t3YXJkKSIseWxhYj0iQWRqdXN0ZWQgUlNxIix0eXBlPSJsIikKcGxvdChiay5zdWJzZXQuc3VtbWFyeSRjcCx4bGFiPSJOdW1iZXIgb2YgVmFyaWFibGVzIChiYWNrd2FyZCkiLHlsYWI9IkNwIix0eXBlPSdsJykKcGxvdChiay5zdWJzZXQuc3VtbWFyeSRiaWMseGxhYj0iTnVtYmVyIG9mIFZhcmlhYmxlcyAoYmFja3dhcmQpIix5bGFiPSJCSUMiLHR5cGU9J2wnKQpgYGAKCmBgYHtyfQpwYXIobWZyb3c9YygyLDIpKQpwbG90KGV4LnN1YnNldC5zdW1tYXJ5JHJzcyx4bGFiPSJOdW1iZXIgb2YgVmFyaWFibGVzIChleGhhdXN0aXZlKSIseWxhYj0iUlNTIix0eXBlPSJsIikKcGxvdChleC5zdWJzZXQuc3VtbWFyeSRhZGpyMix4bGFiPSJOdW1iZXIgb2YgVmFyaWFibGVzIChleGhhdXN0aXZlKSIseWxhYj0iQWRqdXN0ZWQgUlNxIix0eXBlPSJsIikKcGxvdChleC5zdWJzZXQuc3VtbWFyeSRjcCx4bGFiPSJOdW1iZXIgb2YgVmFyaWFibGVzIChleGhhdXN0aXZlKSIseWxhYj0iQ3AiLHR5cGU9J2wnKQpwbG90KGV4LnN1YnNldC5zdW1tYXJ5JGJpYyx4bGFiPSJOdW1iZXIgb2YgVmFyaWFibGVzIChleGhhdXN0aXZlKSIseWxhYj0iQklDIix0eXBlPSdsJykKYGBgCgojIE1vZGVsIEVxdWF0aW9ucwoKSGVyZSB3ZSBwcmludCB0aGUgZmluYWwgZXF1YXRpb25zIGZvciBlYWNoIG1vZGVsLiBOb3QsIHRoZXkgYXJlIGFsbCB0aGUgc2FtZS4KCmBgYHtyfQoKZm9yIChtb2Qub2JqIGluIG1vZGVsLmxpc3QpIHsKICBtb2QubmFtZSA8LSBtb2Qub2JqW1sxXV0KICBiZXN0LmJpYyA8LSBtaW4obW9kLm9ialtbM11dJGJpYykKICBtb2QubnVtIDwtIHdoaWNoLm1pbihtb2Qub2JqW1szXV0kYmljKQogIG1vZC5jYyA8LSBjb2VmKG1vZC5vYmpbWzJdXSwgbW9kLm51bSkKICAKICBtb2QuZXF1YXRpb24uZm9ybWF0IDwtIHBhc3RlKCJZID0iLCBwYXN0ZShyb3VuZChtb2QuY2NbMV0sMiksIAogICAgICAgICAgICAgICAgICAgICAgcGFzdGUocm91bmQobW9kLmNjWy0xXSwyKSwgCiAgICAgICAgICAgICAgICAgICAgICBuYW1lcyhtb2QuY2NbLTFdKSwgCiAgICAgICAgICAgICAgICAgICAgICBzZXA9IiAqICIsIGNvbGxhcHNlPSIgKyAiKSwgCiAgICAgICAgICAgICAgICAgICAgICBzZXA9IiArICIpLCAiKyBlIikKICAKICBwcmludChwYXN0ZSgiTW9kZWwgU2VsZWN0aW9uIE1ldGhvZDogIixtb2QubmFtZSkpCiAgcHJpbnQocGFzdGUoIk1pbiBCSUM6IiwgYmVzdC5iaWMpKQogIHByaW50KHBhc3RlKCJNb2RlbCBOdW1iZXI6ICIsIG1vZC5udW0pKQogIHByaW50KHBhc3RlKCJNb2RlbCBFcXVhdGlvbjogIiwgbW9kLmVxdWF0aW9uLmZvcm1hdCkpCiAgcHJpbnQoIiIpCiAgCn0KYGBgCgoKCgo=